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CLc Abstract 

Based on the weighted and shifted Griinwald difference (WSGD) operators [24], we further construct 
the compact finite difference discretizations for the fractional operators. Then the discretization schemes 
are used to approximate the one and two dimensional space fractional diffusion equations. The detailed 
numerical stability and error analysis are theoretically performed. We theoretically prove and numerically 
1 verify that the provided numerical schemes have the convergent orders 3 in space and 2 in time. 
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1 Introduction 



Diffusion is a fundamental phenomena in real world. The classical diffusion equation can be derived from 



the conservation law of particles, when assuming the particles' diffusion satisfies the Fick's law. Based on 
the CTRW model of statistical physics, it can also be easily derived. Nowadays, anomalous diffusion is 
widely recognized in scientific community, its basic feature is that the classical Fick's law doesn't not hold 
again. In fact, anomalous diffusion is usually characterized by the mean square displacement of the particles, 
given by (x 2 (t)) — (x(t)) 2 ~ t a . The exponent a classifies the type of diffusion: for a = 1, we have normal 
diffusion; for < a < 1, subdiffusion; and a > 1, superdiffusion. The anomalous diffusion equations can 
still be derived from two ways: using the conservation law and fractional Fick's law; and from the CTRW 
model with power law waiting time and/or power law jump length distribution. Fractional calculus plays an 
important role in obtaining the anomalous diffusion equations [1, 2, 3, 4, 10, 11, 19, 25]. 

The equation describing the superdiffusion is the space fractional diffusion equation, and its concrete 
form is to replace the second order derivative of the classical diffusion equation by the Riemann-Liouville 
fractional derivative of order a and 1 < a < 2 [10, 11]. This paper concerns the high accurate numerical 
methods for the space fractional diffusion equation. Meerschaert and his collaborators first introduce the 
shifted Griinwald formula and successfully get the stable finite difference scheme for numerically solving 
the space fractional equation based on this formula [12]. Later more research works are appeared by using the 
shifted Griinwald formula to discretize the space fractional derivatives [8, 12, 13, 14] with first order accuracy 
in space, and possibly obtaining second order accuracy after extrapolation [21, 22]. In [24], the weighted 
and shifted Griinwald difference (WSGD) operators are introduced to discretize the fractional operators 
with second or higher order accuracy, then the second order finite difference schemes are established to 
numerically solve the space fractional diffusion equation. In [24], it is also verified that the 3-WSGD operator 
can approximate the Riemann-Liouville fractional derivative with third order accuracy, but the 3-WSGD 
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operator fails to numerically solve the time dependent space fractional diffusion equations with unconditional 
stability. 

As the sequel to [24], based on the WSGD operators we focus on constructing the compact finite dif- 
ference discretizations, termed compact WSGD operators (CWSGD), for the fractional operators with third 
order accuracy. When the order of fractional derivative a equals to 1 or 2, it becomes the compact difference 
operators for the first or second order spatial derivatives with fourth order accuracy, which have been widely 
used in numerically solving the linear and nonlinear, steady and evolution equations to achieve the high order 
accuracy [7, 17, 20, 23]. By using the CWSGD operator, we get the third order numerical schemes for the 
space fractional diffusion equation. We theoretically make the numerical stability and error analysis, and 
the detailedly numerical experiments confirm the theoretical results. In performing the theoretical analysis, 
the negative definite property of the matrix generated by the WSGD discretization to the fractional diffusion 
operator [24] still plays a key role. 

This paper is organized as follows. In Sec. 2, we introduce the CWSGD operators. Based on the CWSGD 
operators, in Sec. 3 and 4, we build the compact finite difference schemes for the one and two dimensional 
space fractional diffusion equations. The stability and third order convergence with respect to the discrete 
I? norm of the difference schemes are proven. In Sec. 5, the extensive numerical experiments are performed 
to verify the accuracy and theoretical convergent order. And some concluding remarks are made in the last 
section. 

2 Compact WSGD Operators for the Riemami-Liouville Fractional Deriva- 
tives 

Now from the WSGD operator and the Taylor's expansions of the shifted Griinwald finite difference for- 
mulae, we derive the CWSGD operators for the Riemann-Liouville derivatives. First let us introduce the 
definition of the Riemann-Liouville derivatives. 

Definition 1 ([15]). The a(n — 1 < a < n) order left and right Riemann-Liouville fractional derivatives of 
the function u(x) on (a, b) are, respectively, defined as 

(I) left Riemann-Liouville fractional derivative: 



The second order finite difference approximations (the WSGD operators) for Riemann-Liouville frac- 
tional derivatives are given as follows. 

Lemma 1 ([24]). Let u € ^-ooD x +2 u, x D^~ 2 u and their Fourier transformations belong to L 1 (M), 

and define the weighted and shifted Griinwald difference (WSGD) operators by 




(2) right Riemann-Liouville fractional derivative: 




if a = n, then a D%u(x) = ^au(x) and x D£u(x) = (-l) n -^u(x). 



oo 



k=0 
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where g { * ] = {-l) k Q are the coefficients of the power series of the function (1 — z) a , andp, q are integers, 
p 7^ q. Then we have 

L T>l Ptq u(x) = -ooD^x) + 0(h 2 ), (2.2a) 
R Vl Ptq u{x) = x D^u{x) + 0(h 2 ), (2.2b) 

uniformly for x£t 

From [21], we know that the Taylor's expansions of the shifted Griinwald finite difference formulae, 
which is necessary for establishing the CWSGD operator for Riemann-Liouville fractional derivatives. 

Lemma 2 ([21]). Assuming that u £ C n+3 (IR) such that all derivatives of u up to order n + 3 belong to 
then for any nonnegative integer p, we can obtain that 

1 oo n— 1 

— 9^u(x - (k - p)h) = ^D«u(x) + ]T -ooZE+'uOr)) h l + 0(h n ), (2.3a) 
fc=0 l=i 



1 oo n— 1 

^E^ 1 + ( fc " = * D ™< X ) + E (<& xD^ l n(x))h l + 0(h n ), (2.3b) 



fc=0 z=i 
uniformly for x 6 K, where are the coefficients of the power series expansion of function w a)P (z) 



oo 



i=f— ) e*> 2 , and w a Jz) = < k z k = l + (p- f )z + £(a + 3a 2 - 12ap + 12p 2 )z 2 + 0(z 3 ). 

For any function u(x), denoting by J 2 the second order central difference operator, that is S 2 u(x) = 
(u(x — h) — 2u(x) + u(x + h)) /h 2 , we introduce the following finite difference operator 

C x u = (1 + c^ 2 h 2 5 2 x )u, (2.4) 

where c^ q 2 = ^rfyOp, 2 + ^fy«g, 2 - We cal1 C x the CWSGD operator of Riemann-Liouville fractional 
derivatives, of which the detailed construction is described in the following proposition. 

Proposition 1. Under the conditions of Lemmas 1 and 2, there exist 

L Vl M u{x) = C x Li?tt(i)) + t£ ft3 -^D^u^h? + 0(h 4 ), 

) N 7 (2.5) 

rVZ m u(x) = C x ( x D°u(x)) + c£, )3 *r>£+ 3 u(^ 3 + O^ 4 ), 

uniformly for x € R, where p, q are nonnegative integers and p ^ q. The operator C x is defined in (2.4). 

Proof. Substituting formulae (2.3a) and (2.3b) into (2.1a) and (2.1b), respectively, and taking n = 4, we 
arrive at 



L Vl M u{x) = (l + < 3 , 2 /i 2 ^) (-oo^Xx)) + c£, j3 -oo^ +3 «(^)^ 3 + 0(/i 



(2.6) 



Since <S 2 ii = ^it + 0(h 2 ), it yields 



C^ = (l + c^ 2 ^)« + 0(/i 4 ). (2.7) 
Then we obtain the needed results by substituting (2.7) into (2.6). □ 



3 



Remark 1. If the function u(x) is defined on the bounded interval [a, b] with boundary condition u(a) = 
or u(b) = 0, then the WSGD formulae approximating the a order left and right Riemann-Liouville fractional 
derivatives ofu(x) at each point x are written as 



L Vt iPiqU {x) = £ g^ U (x-(k-p)h) + ^ Yl 9^u(x-(k- q )h), 

k=0 k=0 

R Vl m u( X ) = ^ ]T g^ U ( X+ (k- P )h) + ^ Yl 9 { k a) u(x+(k- q )h), 
k=0 k=0 

where )X\ = 2{n-q\ > ^ 2 = 2{p^°q) • Af ter applying Proposition I, we have 
lVI pa u{x) = C x (>>(*)) + c« q>3 a D« +3 u(x)h 3 + 0(/i 4 ), 
RVh, P , q u(x) = C x ( x D?u(x)) + c£ ff>3 x D« +3 u{x)h 3 + 0(h 4 ). 



(2.8) 



(2.9) 



Remark 2. W/je« employing the finite difference method with WSGD formulae for numerically solving non- 
periodic fractional differential equations on bounded interval, p, q should be chosen satisfying \p\ < l,\q\ < 
1 to ensure that the nodes at which the values ofu needed in (2.8) are within the bounded interval; otherwise, 
we need to use another way to discretize the fractional derivative when x is close to the right/left boundary. It 
was indicated in [24] that the approximation by formula (2.8) with (p, q) = (0, —1) turns out to be unstable 
for time dependent problems. Then two sets of (p,q) = (1,0), (1,-1) can be selected to establish the 
difference scheme for fractional diffusion equations, which is also appropriate for the compact difference 
approximations (2.5). The coefficients c pq l in (2.4) with (p, q) = (1,0), (1, —1) are 

Ro,2 = U 7a - 3 « 2 )> c "o,3 = U^~ 3« 2 + 2a), 
K-1,2 = U a ~ 3 « 2 + 12 )> C W,3 = 2i(« 3 " 4a). 

For a = 2, the WSGD operators (2.8) are the centered difference approximation of second order derivative 
when (p,q) equals to (1,0) or (1,-1), and the approximations (2.5) behave as the compact difference 
operators of second order derivative as c\ 2 = c i,-i,2 = 12 an d c i 3 = c i -1 3 = ®> f or a = 1 and 
(p> q) = (1) 0), c\ 2 = g an d c \ 3 = 0» then the centered and compact difference schemes for first order 
derivative are recovered. 

For the cases of (p, q) = (1, 0) and (p, q) = (1, —1), the WSGD schemes at every point X{ = a+i h (h = 
^T, 1 < i < N - 1) are denoted as 



k ~° (2.11) 

N-i+l 



h 

where 



fc=0 



(p >g ) = (1,0), «,W = «,« = f ^ + fc > l; 

(p,g) = (l,-l), w W = 2+£! fl J«), tl ,(-) = 2 + a^), (2.12) 

(a) _ 2 + a (a) 2-a (a) 
w k — ^ 9k + 4 yfc-2' — 
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Remark 3. Let Sn-i be a symmetric tri-diagonal matrix of (N — \)-square, denoted by tridiag(l, —2, 1). 
And we have the eigenvalues of the matrix Sn-i in decreasing order (see [6]), 



A* (5. 



V2iV/ 



N-lj 



k = 1,2,--- ,N-1. 
Define 

C a = In-1 + c p,q,2^N-l, 
where Im-i is the unit matrix of (N — l)-square. Then the eigenvalues ofC a are given by 

kit 



k = 1,2, ••• ,JV-1. 



For the case of (p,q) = (1, 0), we have 

A fe (C Q ) >l-4c? A2 >->0, 

when < a < |; and Xk(C a ) > 1 w/je« a > |. 

Then, from the Rayleigh-Ritz. Theorem (see Theorem 8.8 in [26]), we know that the matrix C D 
c l o 2^N-i) is positive definite. And for the case of(p, q) = (1, —1), we have 

X k {C a ) > 1 - 4 tf ! 2 > iff < a < 



(3 
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(2.13) 
(2.14) 

(2.15) 

(2.16) 
(In-i + 
(2.17) 



1+V73 
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and Afc(C a ) > 1 when a > am/ 1 — 4 c"_ x 2 = w/ze?i a - 

(Jjv— l + c i -i 2^N-i) is positive definite for any natural number N if and only if 'a £ [— 
Lemma 3 ([24]). Le? ?/ze matrix A a be of the following form 



Thus, the matrix C Q 
, oo). 



/ (<*) 



(a) 



(a) 


(a) 
1 

(a) 



\ 



(a) 


(a) 



W 



(a) 
iV-2 
(a) 



.(«) 







(«) ..,(«) 



(2.18) 
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where the diagonals {w k } k =Q are the coefficients given in (2.12) corresponding to (p, q) = (1, 0), (1, —1). 
Then we have that any eigenvalue A of A a satisfies 

(1) Re(X) = 0, for (p,q) = (1,0), a = 1, 

(2) Re(A) < 0,for(p,q) = (1,0), 1 < a < 2, 

(3) Re (A) < 0,for(p,q) = (1,-1), 1 < a < 2. 

Moreover, when 1 < a < 2, matrix A a is negative definite, and the real parts of the eigenvalues A o/ 
matrix c\A a + c 2 ^4q are Zes\s ?/iarc 0, where ci, c 2 > 0, + c 2 7^ 0. 

Lemma 4 ([16, 5]). 77ze Matrix A £ ]R nxn asymptotically stable if and only if there exists a symmetric 
and positive (or negative) definite solution X € W nxn to the Lyapunov equation 

AX + XA T = C, (2.19) 

where C = C T £ M nxn w a negative (or positive) definite matrix. And a matrix A is called asymptotically 
stable if all its eigenvalues have real parts in the open left half -plane, i.e., Re\(A) < 0. 

Lemma 3 and 4 play a central role in analyzing the stability and convergence of the compact difference 
approximations in the sequel. Next taking (p, q) = (1,0) and (p, q) = (1, —1), respectively, we consider the 
compact difference schemes of approximating the space fractional diffusion equations. 
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3 One Dimensional Space Fractional Diffusion Equation 



In this section, we consider the following two-sided one dimensional space fractional diffusion equation 

9 ^f^ = Ki a D«u(x,t) + K 2 x D%u(x,t) + f(x,t), (x,t) e (a,b) x (0,T], 
u(x , 0) = uo(x) , x £ [a,b], ($.1) 

u(a,t) = c(> a (t), u(b,t) = <k(t), t€[0,T], 

where a D% and X D£ are the left and right Riemann-Liouville fractional derivatives with 1 < a < 2, 
respectively. The diffusion coefficients K\ and K 2 are nonnegative constants with Kf + K 2 ^ 0. If K\ ^ 0, 
then 4> a (t) = 0, and if K 2 / 0, then (pb{t) = 0. In the following analysis of the numerical method, we 
assume that (3.1) has a unique and sufficiently smooth solution. 

3.1 Compact Difference Scheme 

We partition the interval [a, b] into a uniform mesh with the space stepsize h = (b — a)/N and the time 
stepsize r = T/M, where N,M are two positive integers. And the set of grid points are denoted by 
Xi = a + ih and t n = nr for < i < N and < n < M. Denoting t n+ i/ 2 = (t n + t n+ \)/2 for 
< n < M — 1, we introduce the following notations 



uf = u(xi,t n ), u™ +1/2 = -(u(xi,t n )+u(xi,t n+ i)), /™ +1 / 2 = f(xi,t n+1/2 ), StUi = « +1 -n™)/r. 
Employing the Crank-Nicolson technique for time discretization of (3.1), we get 

«5 t < - \ (Ki( a D«u)? + K^D^u)^ 1 + K 2 { x D?u)? + K 2 ( X D^ +1 ) = f? +1/2 + 0(r 2 ). (3.2) 



2 

Recalling the definition of operator C x and (2.9), we have 

C x (aD^Ui) = L Vl mUi - c^ 3 h 3 a D^ Ul + 0(/> 4 ), 

v ' (3.3) 
C x ( x D%u>} = RDl >m Ui - c^ 3 h 3 x D%+ 3 Ui + 0(/i 4 ), 

where tD^ v q and rD% are given in (2.11). Acting tC x on both sides (3.2) and then plugging (3.3) into 
it, we obtain 



r n.n , KlT no n , K 2 T n . r ,n+l/2 n+1/2 



(3.4) 



where 



eT 112 = -(K*&fl aD^u^ 1 ' 2 + ^+ 3 < +1/2 > 3 + 0(r* + h% (3.5) 

Denoting by XJf the numerical approximation of u", we derive the compact difference scheme 

„ i+l „ N-i+1 

r U n+1 - — S"w {a) U n+1 V w {a) U n+1 

LxU i Oha Z^ W k U i-k+l ?h a W k U i+k-l 



2h a ^ k 2h c 

k=0 k=0 

-c rr + lT VJ o) fr + 2T V v, {a) rr + T c f n+1/2 

fc=0 fc=0 



(3.6) 
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For the convenience of implementation, we define 



u n = (c/f,c/ 2 v-- M-i) , p n = (/r +1/2 ,/ 2 n+1/2 , 

and reformulate the scheme (3.6) as the following matrix form 



,n+l/2\ T 



C n 



2h c 



■(K 1 A a + K 2 Al))U 



rn+l 



C a + i^{KiA a + K 2 A l a ) )U n + rC a F n + H n , (3.7) 



where y4 a and C a are given in (2.18) and (2.14), respectively, and 



{U n_ U n+l +Tf n^ / , ) + 



n+1/2-. 








"P,Q,2. 



n+l/2x 



2/i c 



K 1 4 a) +Er 2 4 a) ' 



2h c 



K 2 w^ 



K 2 wf ] 
K lW ( a) +K 2 w { 2 a \ 



(3.8) 



3.2 Stability and Convergence 

Next we consider the stability and convergence analysis for the scheme (3.6). Let 

Vh = {v : v = {vi} is a grid function in {xi = a + ih}fL and vo = vn = 0}. 
For any v = {vi} G V/j, we define its pointwise maximum norm and the discrete I? norm as 



w oo = max mj , m 
Ki<AT— 1 



jV-1 
i=l 



(3.9) 



Theorem 1. For the case of (p,q) = (1,0), the difference scheme (3.6) is unconditionally stable for all 
1 < a < 2; and for the case of(p,q) = (1, —1), the difference scheme (3.6) is also unconditionally stable 
for ±±^Zl < a < 2. 

Proo/ Denoting Z) a = i^s(KiA a + K 2 A^), we rewrite (3.7) as 



(C a - D Q )[/ n+1 = (C a + D Q )C/ n + rC Q F n + £T 



n i rrn 



(3.10) 



From Remark 3, we know that C a is a symmetric and positive definite matrix when (p,q) = (1,0) with 

1 < a < 2 and (p, q) = (1,-1) with 1+ ^f^ < a < 2, which follows that C" 1 is also symmetric 

and positive definite. On the other hand, Lemma 3 shows that the eigenvalues of the matrix Da + D °< = 

Ti " K \h^ + ^a) ^ a ^ ne g at i ve for 1 < a < 2, thus (D a + D^) is a symmetric and negative definite 
matrix. Then, for any v = (vi, v 2 , ■ ■ ■ , vn-i) T G M^ -1 \ 0, there exists 



v T ((C^ 1 D a )C~ 1 + C-^C-^afV = v T C-\D a + Dl)C^v < 0, 



(3.11) 
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which means that the matrix ((C Q l D a )C a 1 + C a 1 (C a l D a ) T ) is negative definite. Then it yields from 
Lemma 4 that all the eigenvalues of (C~ l D a ) have negative real parts. In addition, A is an eigenvalue of 
(C" 1 ^) if and only if £±A i s an eigenvalue of matrix (I - C~ l D a )- l {I + C~ l D a ), and |^| < 1 holds. 
Hence, the spectral radius of the matrix (C a — D a )~ l (C a + D a ) = (I — C~ 1 D a )~ 1 (I + C~ 1 D a ) is less 
than one, and the difference scheme (3.6) is stable. □ 

Lemma 5 (Discrete Gronwall Lemma [18]). Assume that {k n } and {p n } are nonnegative sequences, and 
the sequence {4> n } satisfies 

n—l n— 1 

0o < 9o, 4>n < 9o + ^Pi + ^2 h<t>h n>l, 

1=0 1=0 
where qq > 0. Then the sequence {(j> n } satisfies 
n—l n—l 

<f>n < (go + ^2 Pi) ex P {^2 k i\ n>l. (3.12) 

1=0 1=0 

Theorem 2. Let uf be the exact solution of problem (3.1), and Uf be the solution of difference scheme (3.6) 
at grid point (xi,t n ). Then the following estimate 

\\u n - U n \\ < c(t 2 + h 3 ), 1 < n < M, (3.13) 



(3.15) 



holds for all 1 < a < 2 with (p, q) = (1,0) and < a < 2 with (p, q) = (1, -1). 

P roo/ Denoting ef = uf — TJf, from formulae (3.4) and (3.6) we have 

C a (e n+l - e n ) - ^Me n+1 + e n ) - f^I(e n+1 + e n ) = re^\ (3.14) 

where 

e n _ / n n n \T _n+l/2 _ fjn+1/2 n+1/2 n+l/2sT 

e — l e l i e 2 1 " ' > e N-l) ! £ ~~ l e l i e 2 ) J £ JV-1 / • 

Multiplying (3.14) by h(e n+l + e n ) T , we obtain that 

h{e n+1 + e n ) T C a (e n+1 - e n ) - ^rr(e n+1 + e n ) T A a (e n+1 + e n ) 
_ ^^(e^ 1 + e") T ^( e -+i + e ») = r/l ( e n+i + e ")T e n+i/2. 

By Lemma 3, A Q and its transpose are both negative semi-definite matrices for 1 < a < 2, thus 

(e n+1 + e n ) T A a (e n+1 + e n ) < 0, (e n+1 + e n ) T Al(e n+1 + e n ) < 0. (3.16) 
Then (3.15) leads to 

h(e n+1 + e n ) T C a (e n+1 - e n ) < rh(e n+1 + e n ) T e n+1/2 . (3.17) 
As the matrix C a is symmetric, we derive that 

h(e n+1 + e n ) T C a (e n+1 - e n ) = E n+1 - E n , (3.18) 

where 

E n = h(e n ) T C a (e n ) > (1 - 4 c^ 2 )\\e n \\ 2 . (3.19) 
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From (2.16) and (2.17), it yields E n > \\\e n \\ 2 , where A = 1 - 4c?_ 1)2 > if ±±j^Zl < a < 2 and 
(p, = (1, -1); and A = f§ if 1 < a < 2 and (p, q) = (1, 0). Together with (3.17), it yields that 



_ £fc < r/l(e A:+l + e fc)T £ fc+l/2 < I^(|| e fc+l||2 + || e fc||2) + I|| e fc+V2||2. (3 20) 

2 A 

Summing up for all < k < n — 1, we have 

, n— 2 n—2 

X\\e n f < rh(e n + e"" 1 )^"" 1 / 2 + ^ E (H e * +1 || 2 + W^f) + 7 E Ik^ll 2 

fc=0 fc=l 
. 2 n—1 n— 1 

< ^|| e nf + L.|| £ "-V2||2 + rA J- [| e fc||2 + Z_J2 \\e k+1 / 2 \\ 2 . 

k=l k=l 



(3.21) 



to 



(3.22) 



Since |e^ +1//2 | < c(r 2 + h 3 ) for any < k < n — 1, then it leads 
||e"|| 2 < 2r^ ||e fc || 2 + J E H^'ll 2 + ^Ik""^!! 

k=l k=l 
n-l 

<2r^||e fc || 2 +c(r 2 + / l 3 ) 2 , 
k=l 

which completes the proof by Lemma 5. □ 

Remark 4. The truncation error in (3.5) becomes e™ +1/2 = 0(t 2 + h 4 ) when a = 1,2 with (p, g) = (1, 0) 
amc? a = 2 w/?/j (p, g) = (1, —1), so when taking a = 1, 2 f/ie compact finite difference schemes for the 
classical diffusion equations are recovered and the corresponding error estimate of the difference scheme 
(3.6) satisfies 

\\u n -U n \\ <c(r 2 + /i 4 ), l<n<M. (3.23) 

4 Two Dimensional Space Fractional Diffusion Equation 

Now we consider the following two-sided space fractional diffusion equation in two dimensions 

( du{x^t) = (K+ aD>{ ^ ^ t) + K + xD>{Xi Vj t f 

+ (K^ c D^u{x,y,t)^K^ y D p d u{x,y,t)) + f(x,y,t), (x,y,t) G n x (0,T], (41) 

u{x, y, 0) = uo(rc, y), (x, y) G O, 

u(x,y,t) = ip{x,y,t), (x,y,t) G dQ x [0,T], 

where f2 = (a, 6) x (c, d), a D°, x D^ and c Dy, y D^ are Riemann-Liouville fractional derivatives with 
1 < a,f3 < 2. The diffusion coefficients Kf, K~ > 0, i = 1,2, satisfy (i^+) 2 + (K 2 + ) 2 / and 
(Kf) 2 + (i^ - ) 2 / 0. And <p(a, y,t) = if K+ / 0; y>(6, = if / 0; ip(x, c, t) = if K 2 + / 0; 
<p(x, d,t) = if i^2~ 7^ 0- We assume that (4.1) has a unique and sufficiently smooth solution. 
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4.1 Compact Difference Scheme 

We partition the domain Q into a uniform mesh with the space stepsizes h x = (b — a) /N x , h y = (d — c) /N y 
and the time stepsize r = T/M, where N x , N y , M being positive integers. And the set of grid points are 
denoted by X{ = a + ih x ,yj = c + jh y and t n = nr for < i < N x , < j < N y and < n < M. Setting 
tn+i/2 = {tn + t n +i) /2 for < n < M — 1, we denote 

= u{xi,yj,t n ), u™j 1/2 = -(u{xi,yj,t n ) +u{xi,yj,t n+1 )), 

Discretizing (4.1) by the Crank-Nicolson technique in time direction leads to 

SKj = \(Kt{ a D a x u)^ l +K+{ x Dtu)^+K^{ c D^ 



+ K+(aD2u)? d + K+{ x D?u)? d + Jq-( c l6*)?j + ^{yD&Zj ) + /" +1/2 + 0(r 2 ). 



(4.2) 

In the space discretizations, we introduce the finite difference operators 

C x u id = (1 + c^ 2 h 2 x 5l)u id , CyU itj = (1 + 4,g,2 h l S l) u i,3^ ( 4 - 3 ) 
and deduce that 

C K („I>««ij) = /.'P;; - <£ ffj3 /£ tt D«+ 3 « ij - + O(^), (4.4a) 

(xD?^) = flP^^Ui j - < 9)3 ^ x B? +3 Uij + 0{h 4 x ), (4.4b) 

(./>;"-,) = /.p/;,,,,/',, - < 9 ,3^ ^tiy + o(^), (4.40 

= - < g>3 ^ yl^uy + O(^). (4-4d) 



For the simplification of presentation, the same stepsizes are chosen in the following discussions, and denoted 
as h = h x = h y . Acting tC x C v on both sides of (4.2) and plugging (4.4a)-(4.4d) in it, we have 



C * C y h—CyL^h.v.o %— c vR D h.v.a h- C xLT>(_ na ^— C xR V 



.71+1 



2 ~yL u h,p,q 2 V R h,p,q ^ xL h,p,q 2 xR h,p,q) U i,3 
+ TC C f n+1/2 + T£ n+V2 

where 

_n+l/2 _ ,3/^+ „a na+3„, , z^+« r>a+3„, , z^-J? n/3+3 



n+1/2 _ _ , 3 / K + a n a+3 i,+ K + r a D Q+3 7, 4- K~ n/ 3 + ; 



It 



^^4, a , 3y D^ U Y +1/ \o^ + h % 



L 2 ^p,q,3y-^d 

And we denote 

Using the Taylor's expansion shows that 

^^K; 1 -^) = £((J^l£+Jitf,I£)(J^^ (4.7) 
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Adding (4.7) to the corresponding sides of (4.5), then we can factorize it as 

{C ~ \%) {Cy ~ = [C x + T -5i) (Cy + + tCAJ;.; ' 2 + re^ l/ \ (4.8) 



where 

C +V2 = ^^j((«+MM+^)S^+0(^+^ + 4 (4-9) 
By denoting £/Jj as the numerical approximation to ufj, the compact finite difference scheme for (4.1) is 

C x - T -5t) (Cy - ~8l) U^ 1 = (c x + T -5°) (Cy + T -5?) U% + rC x C y f^ 1/2 . (4.10) 

Introducing the intermediate variable V^j, we derive several splitting schemes, such as the compact Peaceman- 
Rachford ADI scheme: 

T \ f T 



C* - 2^JKj = [Cy + -tf)U% + -CyfZ \ (4.11a) 
Cy ~ \ 5 i) U i,r = [C x + \S*)V% + \Cyf-f\ (4.11b) 



the compact Douglas ADI scheme: 

C x - \?>i) V% = (c x C y + r -C y 8« + rC x 8l) U% + rC x C y f^ l/ \ (4. 12a) 

Cy ~ ls£) Ulf 1 = m - l<fe (4.12b) 



and the compact D'Yakonov ADI scheme: 

C x - \&£)V% = (C x + ^) (C, + l*g)t^ + rC,C y /;;/ 1 2 . (4.13a) 



°y- 2 d v) U & =v %- (4 " 13b) 
A simple calculation shows that 

X« W 178 = T ( *^ + K t* D ?)( K I* DP v + K 2yD P d )f^ 1/2 + 0(r 3 h 2 ). (4.14) 
Then from (4.8) and (4.14), it yields that 



2 " X j y-y 2 » y ; «*i j - ^ -r 2 */ V s 2 "/ i,J 

, , r r f n+l/2 r3 *n+l/2 , -n+1/2 

+ TL, x LyJ id + A x°yJi,j +T£ i,j ■ 



(4.15) 



where 



<r /2 = - + K^X^T^ + K 2V^ d )fr; 1/2 + 0(r 2 + r 2 /. 2 ). (4.16) 

Eliminating the truncating error and denoting L/Jj as the numerical approximation of ufj, we have 

^) (c« - p?) t^ 1 = (c, + (c„ + ^ + rc:,c,/-;;; 1 2 + ^ f ;': 1 2 . 

(4.17) 

Introducing the intermediate variable V^ , we obtain the compact locally one-dimensional (LOD) scheme, 

r_\._ r_\._ r/_ r \ n+ i/ 2 



^ - ^s)^ 1 = (c, + T ^] Kj + T -( Cy - T -/ y )^r 12 - 
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4.2 Stability and Convergence 

Let 

Vh = {v : v = {vij} is a grid function in and Vij = on T^}, 

where 

n h = : 1 < i < N x - 1, 1 < j < N y - 1}, 

T h = {(t, j) : i = 0, N x ; < j < N y } U : < i < N x ;j = 0, N y }. 

For any v = {vi} G V^, we define its pointwise maximum norm and discrete L? norm as 



v 



oo 



N x -1 Ny-1 

i=i j=i 



In the following, we first list some properties of Kronecker products of matrices. 

Lemma 6 ([5]). Let A G M nxn /iav<? eigenvalues {Xi}f =v and B G R mxm /j aV e eigenvalues {^j}™^. Then 
the ran eigenvalues of A® B, which represents the kronecker product of matrix A and B, are 

Lemma 7 ([5]). ^ G M mxn , 5 G R rxs , C G R nxp , D G M sx< . Then 

(A ®B)(C®D) = AC® BD (g R mrxpt ). (4.20) 
Moreover, if A, B G ]R™ xri , / a unit matrix of order n, then matrices I A and B ® I commute. 

Lemma 8 ([5]). For all A and B, (A ® B) T = A T ® B T and (A ® B)' 1 = A' 1 ® z/A and 5 are 
/weft/We. 

Lemma 9 ([26]). Le£ ^4, B be two symmetric and positive semi-definite matrices, symbolized A > and 
B>0. Then A®B>0. 

Theorem 3. For the case of (p, q) = (1,0), the difference schemes (4.10) and (4.17) are unconditionally 
stable for 1 < a, j3 < 2. And for the case of (p, q) = (1, —1), the difference schemes (4.10) and (4.17) are 
also unconditionally stable when 1+ ^f^ < ct, (3 < 2. 

Proof. We express grid function C7"- in the vector form as 

Tjn _ l,.n n n n n n n n n \T 

u — V n l,l) "2,1 5 j ^7^-1,1! "1,2) ^2,2) i u N x -l,2i j n l,7V y -l; n 2,7V y -li i u N x -l,N y -l) > 

and denote 



= ® C Q , C y = Cp®I x , (4.21) 



where the symbol ® denotes the Kronecker product, I x and I y are unit matrices of (N x — 1) and (iVy — 1) 
squares, respectively, and the matrices A a and Ap are defined in (2.18) corresponding to a and /3, C a and 
Op are given in (2.14) with coefficients a and /3. Therefore, denoting the disturbances of U n+1 and U n by 
SU n+1 and 5£7™, respectively, we have from (4.10) and (4.17) that 

5U n+1 = (C y -V y y\c x -V x y\c x +V x )(C y + V y )5U n . (4.23) 
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Using Lemma 7, we can check that C x and V x commute with C y and V y , which deduces that (C y — V y ) 1 
and (C y + V y ) commute with (C x — V x ) and (C x + T> x ). Then it obtains from (4.23) that 

5U n = ({C y -V y y\c y + V y )) n ((C x -V x )-\c x +V x )) n 5U°. (4.24) 



From Remark 3, Lemma 6 and Lemma 7, we have that C x and C y are symmetric and positive definite 

matrices in the cases of (p, q) = (1,0) with 1 < a, (3 < 2 and (p, q) = (1,-1) with 1+ ^ < a,f3 < 2, 
which yields that C~ l and C" 1 are also symmetric and positive definite. On the other hand, Lemma 3 and 

8 indicate that the eigenvalues of Ao,J ^ Aa and p+ 2 13 are all negative when 1 < a, /3 < 2, then employing 
Lemma 6, we obtain that (V x +£>J) and (T> y +T> J) are both symmetric and negative definite matrices. Then 
it yields that v T (V x + Pj)t> < and v T (V y + V^)v < hold for any non-zero vector v G ]R( Ar -- 1 )( Ar ^ 1 ), 
and 



T 

V 



(C-^C- 1 + C^iC^V^Av = v T C~ 1 (V 1 + P^CC x t; < 0, j = x,y, (4.25) 



which means that the matrix (CL 1 'D 7 )C 7 1 + C 7 1 (C 7 1 P 7 )' r for 7 = x,y are symmetric and negative 
definite matrices, then it implies from Lemma 4 that the real parts of all the eigenvalues {A 7 } of C" 1 ©-^ for 

7 = x, y are negative, and | jry 1 ! < 1- Additionally, A 7 is an eigenvalue of C^V^ if and only if is 
an eigenvalue of (/ — C~ 1 P 7 )~ 1 (I + C" 1 ^), thus the spectral radius of each matrix is less than 1, which 
concludes that ((/ - C x x V x y l (l + C x x V x )) n and ((J - C y l V y )- l (I + C y l V y )) n converge to zero 
matrix, therefore, the difference scheme (4.10) and (4.17) are stable. □ 

Lemma 10 ([26]). Let A be an n-square symmetric and positive semi-definite matrix. Then there is a unique 
n-square symmetric and positive semi-definite matrix B such that B 2 = A. Such a matrix B is called the 
square root of A, denoted by A 1 / 2 . 

Theorem 4. Let u™j be the exact solution of (4. 1), and f/" - be the solution of the difference schemes (4. 10) or 

(4. 17), then in the cases of (p, q) = (1, 0) with 1 < a, (3 < 2 and (p, q) = (1,-1) with 1+ ^ < a, f3 < 2, 
we have 

\\u n - U n \\ < c(t 2 + h 3 ), 1 < n < M, (4.26) 
where c denotes a positive constant and \\ ■ \\ stands for the discrete L 2 -norm. 
Proof Let efj = uf d - U^, subtracting (4.8) from (4.10) leads to 

(C x - V x ) (C y - Vy)e n+l = (C x + V x ) (C y + V y ) e" + r£ n+1 / 2 , (4.27) 

where 

e = (ei,i, &2,u " " " j e A^-i,ii e i,2, e2,2, • • • ,eN x -i,2,"' > e i,N y -i, e2,N y -i, • • • , e7v ;c _i i 7v y _i) T , 

^ = (^1,1) ^2,1) - - - ,£#3,-1,1) ^1,2, ^2,2) •• • ,£#3,-1,2, •• • ,£l,# y -l, £2,^-1, • • • , £N x -l,N y -l) T , 

and the matrices C x , C y and V x , V y are given by (4.21) and (4.22), respectively. 

As stated in Theorem 3, under the cases of (p, q) = (1, 0) with 1 < a, j3 < 2 and (p, q) = (1, —1) with 
1+ 6^ < ot, P <2, the matrices C Q and Cp and their inverse are symmetric and positive definite. And from 
Lemma 8 and 10, we know that (C" 1 ) 1 / 2 = J v ® (C" 1 ) 1 / 2 and (C" 1 ) 1 / 2 = (C^" 1 ) 1 / 2 <g) I x uniquely exist 
and are symmetric and positive semi-defmite matrices. Then multiplying (4.27) by (C 2 T 1 ) 1 / 2 (C~ 1 ) 1 / 2 , and 
making the discrete L 2 -norm on both sides, we have 

\\(C x 1 ) 1 / 2 (Cy 1 ) 1 / 2 (C x -V x )(C y -Vy)e n+1 \\ 

< \\(C x l ) 1 l\C- 1 ) 1 l\C x + V x )(Cy + Vy)e n \\+T\\(C x ^ 
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Simple calculations show that (C y — V y ) commutes with (C x — T> x ), (C x 1 ) 1 ' 2 , (C x — Pj); (C y + V y ) 
commutes with (C x + V x ), (C" 1 ) 1 / 2 , (C x + Pj); and (Cy 1 ) 1 / 2 commutes with (C" 1 ) 1 / 2 , (C x ± Pj). 
From Lemma 3 and 9, we have that P 7 + P^ (7 = x, y) are symmetric and negative definite matrices. Thus, 
using Lemma 7 and Lemma 9, we obtain that 

((C- 1 ) 1 / 2 ^- 1 ) 1 /^ " 2>.)(C» " -D y )) T ((C- 1 ) 1 / 2 (C- 1 ) 1 /\C x - V x )(C y - V y )) 

= [C y - Vl -V y + V^C^V y ) (C x - pj - P, + PjC"^) ( 4 - 29 ) 
> (C v + -DjC^Vy) (C x + PjC" 1 ^) + (2?J + 2>„) (Pj + P.) , 

and 

((Cx 1 ) 1/2 (C y - 1 ) 1/2 (C7 a! +2?,)(C7 W +2? tf )) T ((C7- 1 ) 1 /2 (C 7-i ) i/2 (C7as + v x )(C y +V v ) j 

= (c y + pj + v y + v^c-'Vy) (c x + v1 + v x + pjc- 1 ©*) ( 4 - 3 °) 
< (a, + vfc-^Vy) (c x + vie- 1 ®*) + (pj + p„) (pJ + v x ) , 

where the matrices A > B means that (A — B) is positive semi-definite. And define 



E n = ^ 2 (e»)T((C„ + VjjC^Vy) (C x + V^C^V X ) + (Pj + p y ) (pT + ^) j ( e n }) (4 . 31) 

it concludes from (4.21), (4.22), Lemma 7, Lemma 9, and Lemma 3 that the matrices C y T>£C~ 1 'D x , 
C x V'£Cy 1 "Dy, PjC-^yPjC- 1 ^ and (Pj + P y ) (Pj + P z ) are all symmetric and positive definite, 
which follows that 



E n > ^^(efCCjfe") = ^h 2 (e™) T (C p ® C a )(e") > A /A mill (C Q )A II1 i I1 (C^)|| e ' l ||, (4.32) 

where A m i n (C a ) and A m i n (C^) are the minimum eigenvalue of matrix C a and Cp, respectively. As stated 
in Remark 3, A min (C a ) > 1 - 4c?_ 1>2 > 0, A min (C /3 ) > 1 - 4cf _ 1>2 > if < a,/3 < 2 and 

(p, (?) = (1, -1); A min (C a ), A min (C^) > f§ if 1 < a, < 2 and (p, g) = (1, 0). Together with (4.29) and 
(4.30), then (4.28) becomes as 

E k+1 _ E k < r || (£-1)1/2^-1)1/2^+1/2 1| _ (4 _ 33) 

From the Rayleigh-Ritz Theorem (see Theorem 8.8 in [26]) and Lemma 6, we have for k = 0, ■ ■ ■ , n — 1 
that 



||(<7-Y /2 (^~Y /2 £ fc+1/2 |l = sjh\Z^ 1 liy(C x 1 C y 1 )E^ 1 l' 1 

< ^KUC^W^W = 1 J |g fc+1/2 1| (4 34) 

£fc+l/2n 



•\/ A m in(Ca)A m i n (C^) 

Summing up (4.33) for all < < n — 1 shows that 

re— 1 Ti—l 

fc=0 V ^min(^aMmin(L-' / 3) 

Combining (4.35) and (4.32), and noticing |e 4 fc j" 1/2 | < c(r 2 + /i 3 ) for all 1 < i,j < N - 1, we obtain 

c 

Amin (Cq ) A m in (C/3 ) 

The estimate for scheme (4. 17) can also be obtained by the similar approach used above. □ 
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Remark 5. If a, (3 = 1, 2 with (p, q) = (1, 0) and a, (3 = 2 vw'f/i (p, q) = (1, —1), f/ien by the reasoning of 
the proof of Theorem 4, we obtain the following error estimate for the difference scheme (4.10) and (4.17) 

\\u n -U n \\ <c(t 2 + h 4 ), l<n<M. (4.37) 

5 Numerical Experiments 

In this section, the numerical results of one and two dimensional cases are presented to show the effectiveness 
and convergence orders of the schemes. 

For saving computational time, we do one extrapolation to increase the accuracy to the third order in 
time (see [9]). The detailed extrapolation algorithm is described as follows. 

Step 1. Calculate £i, (2 from the following linear algebraic equations, 

Ci + C2 = 1, 

< 1 + f = 0; 

Step 2. Compute the solution U n of the compact difference schemes with two time stepsizes r and r/2; 
Step 3. Evaluate the extrapolation solution W u (t) by 

W n (r) = C 1 U n (T) + C 2 U n (r/2). 

Example 1. Consider the following one dimensional space fractional diffusion equation 

= D2u(x, t) + x D?u(x, t) + f(x, t), (x, t) G (0, 1) x (0, 1], 

u(o,t) = «(M) = 0, *e[o,i], (5 - !) 

u{x,0) =x 3 (l-xf, xe [0,1], 



with the source term 



f(x,t) = -e-*(* 3 (l - x) 3 + f ^§- ) (. 3 - + (1 - x) 3 -) - S^^y^ 4 -" + (1 - xf~«) 

+3 W^c7) {x5 ~° + (1 " xf ~ a) " W^) {x6 ~ a + xf ~ a 

And the exact solution of (5.1) is given by u(x, t) = e~ t x 3 (l — x) 3 . 

In Table 1, we present the errors \\u n — W n \\, \\u n — W n \\ OQ and corresponding convergence orders with 
different space stepsizes, where Wp(r) = —\U™{t) + |f/"(r/2) is the extrapolation solution, and C/J 1 
satisfies the compact scheme (3.6). It can be noted that for the case of (p, q) = (1, —1), the numerical results 
are neither stable nor convergent when the order a is less than the critical value 1+ ^ 3 (~ 1.59), which 
coincides with the theoretical results. Figure 1 shows that the convergence rates of the maximum and L 2 
errors to (5. 1) approximated by the compact difference scheme at t = 1 with N = 128 for different a, where 
the convergence rates fall from 4 to 3 near a = 1 and increase gradually with a from 1.9 to 2. 

Example 2. The following fractional diffusion problem 

d^d^ = oD°u(x, y, t) + x D?u(x, y, t) + D^u(x, y, t) + y D^u(x, y, t) + f(x, y, t) (5.2) 

is considered in the domain Q = (0, l) 2 and t > with the boundary conditions 
u(x,y,t) = 0, (x,y)edU, ie[0,l], 



15 



Table 1: The maximum and L 2 errors and corresponding convergence rates to (5.1) approximated by the 
compact difference scheme at t = 1 for different a with r = h. 









(r>, a) = 


(1,0) 






fp, a) = 


(1,-1) 




a 


N 




rate 


\\u n - W n \\ 


rate 


K-Wlloo 


rate 


|| u n -W n \\ 


rate 


1.2 


8 
16 
32 
64 
128 
256 


7.24249E-05 
9.86726E-06 
1.41964E-06 
1.91577E-07 
2.52254E-08 
3.26935E-09 


2.88 
2.80 
2.89 
2.92 
2.95 


4.12739E-05 
6.0055 1E-06 
8.38665E-07 
1.13698E-07 
1.50548E-08 
1.95986E-09 


2.78 
2.84 
2.88 
2.92 
2.94 


1.67141E-01 
6.57473E-01 
1.20582E+04 
6.67837E+11 
6.55394E+25 
7.87651E+50 


-1.98 
-14.16 
-25.72 
-46.48 
-83.31 


1.17292E-01 
4.14981E-01 
7.00900E+03 
3.06749E+11 
2.42775E+25 
2.34722E+50 


-1.82 
-14.04 
-25.38 
-46.17 
-83.00 


1+V73 
6 


8 

16 

32 
64 
128 
256 


5.48070E-05 
6.71566E-06 
8.92036E-07 
1.19966E-07 
1.60292E-08 
2.11675E-09 


3.03 
2.91 
2.89 
2.90 
2.92 


3.18317E-05 
4.15788E-06 
5.69588E-07 
7.81948E-08 
1.06194E-08 
1.42376E-09 


2.94 
2.87 
2.86 
2.88 
2.90 


2.66168E-04 
4.38053E-05 
6.07692E-06 
8.00832E-07 
1.05299E-07 
1.37085E-08 


2.60 
2.85 
2.92 
2.93 
2.94 


1.62093E-04 
2.50890E-05 
3.69218E-06 
5.16694E-07 
7.0025 1E-08 
9.30500E-09 


2.69 
2.76 
2.84 
2.88 
2.91 



1.8 



8 


3.78749E-05 




2.87085E-05 




1.49798E-04 




8.54691E-05 




16 


3.90569E-06 


3.28 


2.91599E-06 


3.30 


1.94926E-05 


2.94 


1.23638E-05 


2.79 


32 


4.51151E-07 


3.11 


3.47233E-07 


3.07 


2.60919E-06 


2.90 


1.76317E-06 


2.81 


64 


5.86615E-08 


2.94 


4.44726E-08 


2.96 


3.31182E-07 


2.98 


2.42232E-07 


2.86 


128 


7.68726E-09 


2.93 


5.85045E-09 


2.93 


4.26392E-08 


2.96 


3.24700E-08 


2.90 


256 


1.00982E-09 


2.93 


7.74722E-10 


2.92 


5.51021E-09 


2.95 


4.289 18E-09 


2.92 




1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 '1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 



a a 
(a) (p,<z) = (l,0) (b) (p,q) = (1,-1) 

Figure 1: The convergence rates of the maximum and I? errors to (5.1) approximated by the compact 
difference scheme at t = 1 with N = 128. 

and initial value 

u(x, y, 0) = u(x, y, 0) = x 3 (l - x) V (1 - yf, (x, y) G [0, l] 2 . 
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The source term is 



f{x,y,t) 



x 3 (l - x) 3 y 3 (l - y) 



+ 
+ 
+ 
+ 



r(4) 



„3-a 



„5— a 



T(4 - a) 

3r(6) , 

r(6 - a) 1 

r(4) 
r(4-/3) 

3r (6) t..s-p 



+ (l-xf- a ) 



+ (1 -x) 5 ^ Q ) 



^ + (1 - yf-V) 



3r(5) 

r(5-a 

r(7) 

T(7-a) 
3T(5) 



(x 4 " a + (1 - x) 4 " a ) 



r(5-/3) 



(a; 6 - a + (l-x) 6 ^))y 3 (l-y) J 

(y 4 ^ + (l - y) 4 "' 3 ) 



y^P + (1 _ y) 5^) 



r(7) 



r(6-/?) vy ■ v- " / r(7-/3) 

A«(i exacf solution of (5.1) is given by u(x,t) = e~ t x s (l — x) 3 y 3 (l — y) 3 . 



( y 6-^ + ( l_ y) 6-^ U3 (1 _ x) 



Table 2: The maximum and L 2 errors and corresponding convergence rates to (5.1) approximated by the 
compact difference splitting schemes at t = 1 with r = h. 







(p. n) — 


(1,0), (a 


,8) = (1.1, 1.7) 




(»,o) = (1, 




a,0) = (1.6,1.9) 




Scheme 


A T 

A 


\\vr - V^™||oo 


rate 


1 77 T J 7 71, 1 1 

\\u — W || 


rate 


\\U n - W n \\oo 


rate 


1 77 TI 7"77 1 1 

| \u — W || 


rate 




8 


5.18337E-07 




2.35962E-07 




1.27051E-05 




5.20284E-06 






16 


7.47852E-08 


2.79 


3.12736E-08 


2.92 


8.45409E-07 


3.91 


3.09460E-07 


4.07 




32 


9.61857E-09 


2.96 


4.13976E-09 


2.92 


8.71092E-08 


3.28 


2.87346E-008 


3.43 


CLOD 


64 


1.28508E-09 


2.90 


5.56504E-10 


2.90 


9.67026E-09 


3.17 


3.35453E-009 


3.10 




128 


1.71211E-10 


2.91 


7.47594E-11 


2.90 


1.17863E-09 


3.04 


4.22753E-010 


2.99 




256 


2.26863E-11 


2.92 


9.97266E-12 


2.91 


1.49079E-10 


2.98 


5.43715E-011 


2.96 




8 


6.30467E-007 




2.63036E-007 




3.59941E-006 




1.33801E-006 






16 


7.67884E-008 


3.04 


3.20134E-008 


3.04 


5.12851E-007 


2.81 


1.88287E-007 


2.83 




32 


9.54417E-009 


3.01 


4.19273E-009 


2.93 


7.24724E-008 


2.82 


2.48958E-008 


2.92 


CPR 


64 


1.28241E-009 


2.90 


5.60865E-010 


2.90 


9.12654E-009 


2.99 


3.25063E-009 


2.94 




128 


1.71042E-010 


2.91 


7.51181E-011 


2.90 


1.15503E-009 


2.98 


4.22913E-010 


2.94 




256 


2.26759E-011 


2.92 


1.00029E-011 


2.91 


1.49068E-010 


2.95 


5.48530E-011 


2.95 




8 


6.28561E-007 




2.54619E-007 




3.01132E-006 




1.15694E-006 






16 


7.67367E-008 


3.03 


3.12227E-008 


3.03 


4.37299E-007 


2.78 


1.67892E-007 


2.78 




32 


9.54020E-009 


3.01 


4.12678E-009 


2.92 


6.67372E-008 


2.71 


2.30759E-008 


2.86 


CDouglas 


64 


1.28220E-009 


2.90 


5.557 10E-010 


2.89 


8.61113E-009 


2.95 


3.09798E-009 


2.90 




128 


1.71028E-010 


2.91 


7.47144E-011 


2.89 


1.11175E-009 


2.95 


4.09795E-010 


2.92 




256 


2.26748E-011 


2.92 


9.97005E-012 


2.91 


1.45187E-010 


2.94 


5.36563E-011 


2.93 




8 


6.28561E-007 




2.54619E-007 




3.01132E-006 




1.15694E-006 






16 


7.67367E-008 


3.03 


3.12227E-008 


3.03 


4.37299E-007 


2.78 


1.67892E-007 


2.78 




32 


9.54020E-009 


3.01 


4.12678E-009 


2.92 


6.67372E-008 


2.71 


2.30759E-008 


2.86 


CD'yakonov 


64 


1.28220E-009 


2.90 


5.55710E-010 


2.89 


8.61113E-009 


2.95 


3.09798E-009 


2.90 




128 


1.71028E-010 


2.91 


7.47144E-011 


2.89 


1.11175E-009 


2.95 


4.09795E-010 


2.92 




256 


2.26748E-011 


2.92 


9.97005E-012 


2.91 


1.45187E-010 


2.94 


5.36563E-011 


2.93 



In Table 2, the errors \\u n — W n \\, \\u n — TV n ||oo and their respective convergence rates are presented 
for different uniformly space stepsizes, where WPj(r) = —^U^^t) + |f7 l "-(r/2) is the numerical solution 
by extrapolation in time, and [/J 1 - satisfies the compact LOD scheme (4.18), compact Peaceman-Richardson 
scheme (4.11), compact Douglas scheme (4.12) and compact D'yakonov scheme (4.13), respectively. The 
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third order accuracy both in time and space is verified, and in the computational process, the time costs are 
largely reduced. 

6 Conclusion 

In [24], we introduce the weighted and shifted Griinwald difference (WSGD) operators and show that the 
WSGD operators have second order accuracy to approximate the fractional derivatives. This paper is the 
sequel of [24]. Based on the WSGD operators, we further introduce the compact WSGD operators (CWSGD) 
which have third order accuracy. Then we use the CWSGD operators to establish the compact difference 
schemes for one and two dimensional space fractional diffusion equations. And the theoretical analysis of 
the stability and convergence of the schemes is presented. The numerical results illustrate the effectiveness 
of the compact difference approximation for the fractional problems and confirm the convergence orders of 
the schemes. 
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